#include <time.h>
#include <math.h>
#include "eno.h"
#include "fluids.h"

#include "eno.xtn"
#include "fluids.xtn"

extern int number_time_steps;

extern double compute_curvature(),compute_delta(),
              compute_delta_general(),compute_heaviside(),
              compute_heaviside_general(),compute_delta_prime();
extern double b00[imax][jmax],b01[imax][jmax],b02[imax][jmax],
              b10[imax][jmax],b11[imax][jmax],b12[imax][jmax];

double global_area;

static double mass_spread;
static int i_dim;

#define plotmax (1000)

double x1plot[plotmax],y1plot[plotmax],x2plot[plotmax],y2plot[plotmax];
int plot_count;

double xc3(i,iub,xc)
int i,iub;
double xc;
{
  if (i==0)
   return xmin;
  else if (i==iub+1)
   return xmax;
  else return xc;
}

double yc3(j,jub,yc)
int j,jub;
double yc;
{
  if (j==0)
   return ymin;
  else if (j==jub+1)
   return ymax;
  else return yc;
}

#define staggered_grid 1
#define stag ( (staggered_grid) ? 1 : 0)
#define xc1(i) ( xmin+((i)-1)*hx+0.5*hx*stag )
#define yc1(j) ( ymin+((j)-1)*hy+0.5*hy*stag)
#define xc(i) ( xc3((i),iub,xc1(i)) )
#define yc(j) ( yc3((j),jub,yc1(j)) )
#define zz1(i,j) ( gddd(dd,(i),(j)) )

void sum_global_area(x1,y1,x2,y2,dd,iub,jub,hx,hy)
int iub,jub;
double x1,y1,x2,y2,hx,hy;
double dd[imax][jmax];
{
double mm,parm,y2end,y1end,axi_area,xs,ys,hz[2][2],
       xd[2],yd[2],hs[2],havg,xmid;
int ix,iy,ix1,iy1;

  x1plot[plot_count]=x1;x2plot[plot_count]=x2;
  y1plot[plot_count]=y1;y2plot[plot_count]=y2;
  if (plot_count+1<plotmax)
   plot_count++;

  if ((use_symmetric)||(use_axisymmetric))
   xmid=0.0;
  else
   xmid=(xmax+xmin)/2.0;

  if ((fabs(y2-y1)>almost_same)&&(x1>=xmid)&&(x2>=xmid))
    {
     mm=(x2-x1)/(y2-y1);
     parm=x1-mm*y1-xmid;

     if (use_axisymmetric)
      {
       y2end=sqr(parm)*y2+mm*sqr(y2)*parm+sqr(mm)*sqr(y2)*y2/3.0;
       y1end=sqr(parm)*y1+mm*sqr(y1)*parm+sqr(mm)*sqr(y1)*y1/3.0;
       axi_area=pi*fabs(y2end-y1end);
      }
     else
      {
       y2end=parm*y2+mm*sqr(y2)/2.0;
       y1end=parm*y1+mm*sqr(y1)/2.0;
       axi_area=2.0*fabs(y2end-y1end); 
       /* factor of 2 due to assumed symmetry */
      }

     if (fabs(x1-x2)<sqr(hx))
      {
       xs=max(x1,x2)-sqr(hx);
       if (x1<x2) ys=y2;
       else ys=y1;
      }
     else if (x1<x2) 
      { 
       xs=x1;ys=y2; 
      }
     else 
      { 
       xs=x2;ys=y1; 
      }

     ix=(int) ( floor( (xs-xmin-hx/2.0)/hx )+1 );
     iy=(int) ( floor( (ys-ymin-hy/2.0)/hy )+1 );
     for (ix1=0;ix1<2;ix1++)
      for (iy1=0;iy1<2;iy1++)
       hz[ix1][iy1]=zz1(ix1+ix,iy1+iy);

     xd[0]=xc(ix)-xs;xd[1]=xc(ix+1)-xs;
     yd[0]=yc(iy)-ys;yd[1]=yc(iy+1)-ys;

     for (iy1=0;iy1<2;iy1++)
      hs[iy1]=(xd[1]*hz[0][iy1]-xd[0]*hz[1][iy1])/(xd[1]-xd[0]);
     havg=(yd[1]*hs[0]-yd[0]*hs[1])/(yd[1]-yd[0]);

     if (havg>0.0)
      axi_area=-axi_area;

if ((havg>0.0)&&(FALSE))
{
printf("x1,y1,x2,y2,xs,ys%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n",
 x1,y1,x2,y2,xs,ys);
printf("ix,iy,havg,axi_area%5d%5d%15.12f%10.5f\n",ix,iy,havg,axi_area);
}

     global_area += axi_area;
    }    /* no area/volume for horizontal lines */
}    /* sum_global_area */


void zero_level_set(iub,jub,dd)
int iub,jub;
double dd[imax][jmax];
{
double zc=0.0;
double h[5],xh[5],yh[5];
int castab1[27];
int castab[3][3][3];
int ish[5];
int im[4];
int jm[4];
int i,j,k,l,m,m1,m2,m3,cse;
double x1,y1,x2,y2,dmin,dmax,hx,hy;

  plot_count=0;

  im[0]=0;im[1]=1;im[2]=1;im[3]=0;
  jm[0]=0;jm[1]=0;jm[2]=1;jm[3]=1;
  castab1[0]=0;castab1[1]=0;castab1[2]=8;
  castab1[3]=0;castab1[4]=2;castab1[5]=5;
  castab1[6]=7;castab1[7]=6;castab1[8]=9;

  castab1[9]=0;castab1[10]=3;castab1[11]=4;
  castab1[12]=1;castab1[13]=3;castab1[14]=1;
  castab1[15]=4;castab1[16]=3;castab1[17]=0;

  castab1[18]=9;castab1[19]=6;castab1[20]=7;
  castab1[21]=5;castab1[22]=2;castab1[23]=0;
  castab1[24]=8;castab1[25]=0;castab1[26]=0;

  l=0;
  for (k=0;k<3;k++)
   for (j=0;j<3;j++)
    for (i=0;i<3;i++)
     {
      castab[k][j][i]=castab1[l];
      l++;
     } 

  hx=(xmax-xmin)/(iub-1.0+stag);
  hy=(ymax-ymin)/(jub-1.0+stag);

  global_area=0.0;

  for (j=jub-1+stag;j>=1-stag;j--)
   for (i=1-stag;i<iub+stag;i++)
    {

     /* kluge for computing area where fluid goes up against rt wall */
     h[2]=zz1(i+1,j);h[3]=zz1(i+1,j+1);
     xh[2]=xc(i+1);yh[2]=yc(j);
     xh[3]=xc(i+1);yh[3]=yc(j+1);

     if (i==iub)
      if ((h[2]<=zc)&&(h[3]<=zc))
       sum_global_area(xh[2],yh[2],xh[3],yh[3],dd,iub,jub,hx,hy);
      else
       {
        if (fabs(h[2]-h[3])>almost_same)
         y2=(h[2]*yh[3]-h[3]*yh[2])/(h[2]-h[3]);
        else
         y2=yh[2];

        if (h[2]<=zc)
         sum_global_area(xh[2],yh[2],xh[3],y2,dd,iub,jub,hx,hy);
        else if (h[3]<=zc)
         sum_global_area(xh[3],yh[3],xh[2],y2,dd,iub,jub,hx,hy);
       }

     dmin=min(zz1(i,j),zz1(i,j+1));
     dmin=min(dmin,zz1(i+1,j));
     dmin=min(dmin,zz1(i+1,j+1));

     dmax=max(zz1(i,j),zz1(i,j+1));
     dmax=max(dmax,zz1(i+1,j));
     dmax=max(dmax,zz1(i+1,j+1));
    
     if ((zc>=dmin) && (zc<=dmax))
      {
       for (m=4;m>=0;m--)
        {
         if (m>0)
          {
           h[m]=zz1(i+im[m-1],j+jm[m-1])-zc;
           xh[m]=xc(i+im[m-1]);yh[m]=yc(j+jm[m-1]);
          }
         if (m==0)
          {
           h[0]=(h[1]+h[2]+h[3]+h[4])/4.0;
           xh[0]=(xc(i)+xc(i+1))/2.0;yh[0]=(yc(j)+yc(j+1))/2.0;
          }
         if (h[m]>0.0) ish[m]=2;
         else ish[m]=( (h[m]<0.0) ? 0 : 1 );
        }   /* looping m */

 
       for (m=1;m<=4;m++)
        {
         m1=m;m2=0;m3=m+1;
         if (m3==5) m3=1;

         cse=castab[ ish[m1] ][ ish[m2] ][ ish[m3] ];

         if (cse==1)
          {
           x1=xh[m1];
           y1=yh[m1];
           x2=xh[m2];
           y2=yh[m2];
          }
         if (cse==2)
          {
           x1=xh[m2];
           y1=yh[m2];
           x2=xh[m3];
           y2=yh[m3];
          }
         if (cse==3)
          {
           x1=xh[m3];
           y1=yh[m3];
           x2=xh[m1];
           y2=yh[m1];
          }
         if (cse==4)
          {
           x1=xh[m1];
           y1=yh[m1];
           x2=(h[m3]*xh[m2]-h[m2]*xh[m3])/(h[m3]-h[m2]);
           y2=(h[m3]*yh[m2]-h[m2]*yh[m3])/(h[m3]-h[m2]);
          }
         if (cse==5)
          {
           x1=xh[m2];
           y1=yh[m2];
           x2=(h[m1]*xh[m3]-h[m3]*xh[m1])/(h[m1]-h[m3]);
           y2=(h[m1]*yh[m3]-h[m3]*yh[m1])/(h[m1]-h[m3]);
          }
         if (cse==6)
          {
           x1=xh[m3];
           y1=yh[m3];
           x2=(h[m2]*xh[m1]-h[m1]*xh[m2])/(h[m2]-h[m1]);
           y2=(h[m2]*yh[m1]-h[m1]*yh[m2])/(h[m2]-h[m1]);
          }
         if (cse==7)
          {
           x1=(h[m2]*xh[m1]-h[m1]*xh[m2])/(h[m2]-h[m1]);
           y1=(h[m2]*yh[m1]-h[m1]*yh[m2])/(h[m2]-h[m1]);
           x2=(h[m3]*xh[m2]-h[m2]*xh[m3])/(h[m3]-h[m2]);
           y2=(h[m3]*yh[m2]-h[m2]*yh[m3])/(h[m3]-h[m2]);
          }
         if (cse==8)
          {
           x1=(h[m3]*xh[m2]-h[m2]*xh[m3])/(h[m3]-h[m2]);
           y1=(h[m3]*yh[m2]-h[m2]*yh[m3])/(h[m3]-h[m2]);
           x2=(h[m1]*xh[m3]-h[m3]*xh[m1])/(h[m1]-h[m3]);
           y2=(h[m1]*yh[m3]-h[m3]*yh[m1])/(h[m1]-h[m3]);
          }
         if (cse==9)
          {
           x1=(h[m1]*xh[m3]-h[m3]*xh[m1])/(h[m1]-h[m3]);
           y1=(h[m1]*yh[m3]-h[m3]*yh[m1])/(h[m1]-h[m3]);
           x2=(h[m2]*xh[m1]-h[m1]*xh[m2])/(h[m2]-h[m1]);
           y2=(h[m2]*yh[m1]-h[m1]*yh[m2])/(h[m2]-h[m1]);
          }
         if (cse>0)
          {
           sum_global_area(x1,y1,x2,y2,dd,iub,jub,hx,hy);
          }
        }  /* looping m - all 4 triangles */
      }   /* square had contour line through it */
    }   /* looping all squares on the mesh */

}   /* zero_level_set */

/* compute area that the level set function dd is less than zero.  */
void get_current_area(out_area,dd)
double *out_area,dd[imax][jmax];
{
int area_mx,area_my,ix1,ix2,iy1,iy2,ix3,iy3,ix,iy,ix_alt,iy_alt;
double area_h,x1,x2,x3,y1,y2,y3,hld,sum_pos,sum_tot,
  d1,d2,d3,d4,d5;
double var[3][3],r_minus,r_plus,r_val;
clock_t t0,t1;

#define old_area FALSE
#define fast_area FALSE

  t0=clock();

if (old_area)
{
 zero_level_set(mx,my,dd);
 *out_area=global_area;
}
else if (fast_area)
{
 for_all_points
   if (! point_ok)
    dd[ix][iy]=gddd(dd,ix,iy);

 *out_area=0.0;
 r_val=r_plus=r_minus=1.0;
 for_all_points
  if (point_ok)
   {
    if (use_axisymmetric)
     {
      r_val=x[ix];
      r_plus=r_val+h;
      r_minus=r_val-h;  /* we allow r_minus<0 */
     }
    for (ix1=0;ix1<3;ix1++)
     for (iy1=0;iy1<3;iy1++)
      {
       var[ix1][iy1]=0.0;
       if (dd[ix+ix1-1][iy+iy1-1]<=0.0)
         var[ix1][iy1]=1.0;
       }
     *out_area+=( (r_minus*(var[0][0]+var[0][2])+r_plus*(var[2][2]+var[2][0])+
       22*(r_minus*var[0][1]+r_plus*var[2][1]+r_val*(var[1][0]+var[1][2]))+
       484*r_val*var[1][1])*sqr(h)/576.0 );
   }  /* looping all points */

 if (use_axisymmetric)
  *out_area*=(2.0*pi);
 else if (use_symmetric)
  *out_area*=2.0;
}   /* above if (fast_area) */
else
{
#define area_m 10

  area_h=h/area_m;

  *out_area=0;
  for_all_points
   if (! point_ok)
    dd[ix][iy]=gddd(dd,ix,iy);
  x[0]=x[1]-h;y[0]=y[1]-h;
  x[mx+1]=x[mx]+h;y[mx+1]=y[mx]+h;

#define interp(ixval,iyval,xval,yval) \
 (  ( ((x[ixval+1]-xval)/h)*dd[ixval][iyval]+  \
      ((xval-x[ixval])/h)*dd[ixval+1][iyval] )*((y[iyval+1]-yval)/h) +  \
    ( ((x[ixval+1]-xval)/h)*dd[ixval][iyval+1]+  \
      ((xval-x[ixval])/h)*dd[ixval+1][iyval+1] )*((yval-y[iyval])/h)  )

  for_domain
   {
    sum_pos=sum_tot=0.0;
    for (ix_alt=ix-1;ix_alt<=ix+1;ix_alt++)
     for (iy_alt=iy-1;iy_alt<=iy+1;iy_alt++)
      {
       hld=dd[ix_alt][iy_alt];
       if (hld>0.0)
        {
         sum_pos+=hld;
         sum_tot+=hld;
        }
       else
        sum_tot-=hld;
      }
    if (sum_pos==sum_tot)
     hld=0.0;
    else if (sum_pos==0.0)
     {
      if (use_axisymmetric)
       *out_area+=pi*(sqr(x[ix]+h/2)-sqr(x[ix]-h/2))*h;
      else if (use_symmetric)
       *out_area+=sqr(h)*2.0;
      else
       *out_area+=sqr(h);
     }
    else
     {
     for (ix_alt=1;ix_alt<=area_m;ix_alt++)
     for (iy_alt=1;iy_alt<=area_m;iy_alt++)
     {
     x1=x[ix]-h/2+(ix_alt-1)*area_h;
     x2=x1+area_h;x3=x1+area_h/2;
     y1=y[iy]-h/2+(iy_alt-1)*area_h;
     y2=y1+area_h;y3=y1+area_h/2;
     ix1=( (x1<x[ix]) ? ix-1 : ix);
     ix2=( (x2<x[ix]) ? ix-1 : ix);
     ix3=( (x3<x[ix]) ? ix-1 : ix);

     iy1=( (y1<y[iy]) ? iy-1 : iy);
     iy2=( (y2<y[iy]) ? iy-1 : iy);
     iy3=( (y3<y[iy]) ? iy-1 : iy);

     d1=interp(ix1,iy1,x1,y1);
     d2=interp(ix2,iy1,x2,y1);
     d3=interp(ix1,iy2,x1,y2);
     d4=interp(ix2,iy2,x2,y2);
     d5=interp(ix3,iy3,x3,y3);
     sum_pos=(max(d1,0.0)+max(d2,0.0)+max(d3,0.0)+max(d4,0.0)+4*max(d5,0.0));
     sum_tot=fabs(d1)+fabs(d2)+fabs(d3)+fabs(d4)+4*fabs(d5);
     if (sum_tot==0.0)
      hld=0.0;
     else
      hld=(sum_tot-sum_pos)/sum_tot;

     if (use_axisymmetric)
      *out_area+=pi*(sqr(x2)-sqr(x1))*area_h*hld;
     else if (use_symmetric)
      *out_area+=2.0*sqr(area_h)*hld;
     else
      *out_area+=sqr(area_h)*hld;
     }   /* looping inside the rectangle */
     }  /* had to disect the rectangle */
   } /* looping the domain */
}  /* doing new not old */

  t1=clock();
/*  printf("area,time %10.6f %10ld\n",*out_area,t1-t0);*/
     
}

/* resolve cases in which all 3 level set functions 
   are either negative or positive */
void handle_triple(d1in,d2in,d3in,cutoff)
double *d1in,*d2in,*d3in,cutoff;
{
double d1,d2,d3,d1n,d2n,d3n,dmin,dmax,dmid,dhold;

if ((three_phase())&&
    ((feed_back()) || (reinit_all))
   )
 {
  d1 = *d1in;
  if (! feed_back() )
   d1=-d1;
  d2 = *d2in;
  d3 = *d3in;

  d1n = 0.5*(d1-min(d2,d3));
  d2n = 0.5*(d2-min(d1,d3));
  d3n = 0.5*(d3-min(d1,d2));

#define swap_d(dfirst,dsecond) if (dfirst>dsecond) {dhold=dfirst; \
    dfirst=dsecond;dsecond=dhold;}
#define use_bmo_equation 1

  dmin=d1;
  dmid=d2;
  dmax=d3;
  swap_d(dmin,dmid)
  swap_d(dmid,dmax)
  swap_d(dmin,dmid)

  if (feed_back() )
   {
    if (use_bmo_equation)
     {
      *d1in=d1n;
      *d2in=d2n;
      *d3in=d3n;
     }
    else
     {
      if ((d1<=dmid)||(dmax<=0.0))
       *d1in=d1n;
      if ((d2<=dmid)||(dmax<=0.0))
       *d2in=d2n;
      if ((d3<=dmid)||(dmax<=0.0))
       *d3in=d3n;
     }
   }  /* if feed_back() */
  else
   {
    if (d1>dmid+almost_same)
     {
      *d2in=d2n;
      *d3in=d3n;
     }
    else if (d2<=dmid)
     *d2in=-d1;
    else if (d3<=dmid) 
     *d3in=-d1;
   }    /* handling triple point when no feed back */
 }  /* if three_phase() */
}   /* handle_triple */

void resolve_triple_point()
{
int ix,iy;
double cutoff;

 if ((three_phase())&&
     ((feed_back()) || (reinit_all))
    )
  {
   cutoff = density_spread*h;

   for_all_points
    if (point_ok)
      handle_triple(&d1[ix][iy],&d2[ix][iy],&d3[ix][iy],cutoff);
  }   /* if three_phase() */
}  /* resolve_triple_point */

void calibrate()
{
int ix,iy;
double c1,c2,c3,h1,h2,h3;

  for_all_points
   if (point_ok)
    {
     c1=d1[ix][iy];c2=d2[ix][iy];c3=d3[ix][iy];
     if ((c1>=c2)&&(c1>=c3))
      {
       h1=c1;h2=(c2-c3)/2.0;h3=-c1;
      }
     else if ((c2>=c1)&&(c2>=c3))
      {
       h1=-c2;h2=c2;h3=(c3-c1)/2.0;
      }
     else if ((c3>=c1)&&(c3>=c2))
      {
       h1=(c1-c2)/2.0;h2=-c3;h3=c3;
      }

     d1[ix][iy]=h1;
     d2[ix][iy]=h2;
     d3[ix][iy]=h3;
    }
}   /* calibrate */

void uncalibrate()
{
int ix,iy;
double c1,c2,c3;

  for_all_points
   if (point_ok)
    {
     c1=d1[ix][iy];c2=d2[ix][iy];c3=d3[ix][iy];
     d1[ix][iy]=max(c1,-c3);
     d2[ix][iy]=max(c2,-c1);
     d3[ix][iy]=max(c3,-c2);
    }
}   /* uncalibrate */

#define old_dd(ix,iy) (dd[ix][iy])
#define new_dd(ix,iy) (d4[ix][iy])
#define new_dd1(ix,iy) ( gddd(d4,(ix),(iy)) )
#define grad_dd(ix,iy) (yy[ix][iy])
#define old_grad_dd(ix,iy) (d6[ix][iy])
#define time_correct (rk_time_reinit)
#define base_correct(ix,iy) (old_dd(ix,iy))
#define dt_correct(ix,iy) ((new_dd(ix,iy)-base_correct(ix,iy))/time_correct)
#define grad_correct(ix,iy) (old_grad_dd(ix,iy))
#define delta_correct(ix,iy) \
 (compute_delta_general(old_dd(ix,iy),mass_spread))

#define rk_max 2

double compute_bias(main_level,bias_level)
double main_level,*bias_level;
{
  *bias_level=main_level;
  return (2.0*compute_heaviside_general(*bias_level,mass_spread));
/* we use mass_spread here, NOT density_spread, density_spread only
   used for determining the iteration stop time */

}

/* reinitialize the level set function dd using an iterative eno-type scheme.
*/
void check_mass_conserve(dd,rk_time_reinit)
double dd[imax][jmax],rk_time_reinit;
{
int ix,iy,icx,icy;
double new_level,current_level,main_level,grad_temp,new_teemp,dt_temp,
  delta_temp,bl_int[3][3],cl_int[3][3],cl_sum,bl_sum,f_temp,f_center;
double lambda,lambda_linear,r_val,r_plus,r_minus,dchange,d7_hold;


/* correct d4 to conserve mass - temporary correction in d7 */ 

  for_all_points
   {
    if (! point_ok)
     {
      d4[ix][iy]=gddd(d4,ix,iy);
      d6[ix][iy]=gddd(d6,ix,iy);
      yy[ix][iy]=gddd(yy,ix,iy);
     }  
    d7[ix][iy]=d4[ix][iy];
   }

  for_domain
   {
    f_center=delta_correct(ix,iy)*grad_correct(ix,iy); 
    /* delta(0)  NOT h'(0)  */

    if (fabs(f_center)>almost_same)
     {
      current_level = base_correct(ix,iy);  /* p0 */
      dchange=dt_correct(ix,iy);     /*  (pn-p0)/time */

      r_val=r_plus=r_minus=1.0;
      if (use_axisymmetric)
       {
        r_val=x[ix];
        r_plus=r_val+h/4.0;
        r_minus=r_val-h/4.0; 
       }

      for (icx=0;icx<3;icx++)
       for (icy=0;icy<3;icy++)
        {
         delta_temp=delta_correct(ix+icx-1,iy+icy-1); 
         dt_temp=dt_correct(ix+icx-1,iy+icy-1); 
         /* (pn-p0)/time */

         f_temp=delta_temp*grad_correct(ix+icx-1,iy+icy-1);
         /* delta(0)  NOT h'(0)  */

         bl_int[icx][icy]=-f_temp*delta_temp;
         cl_int[icx][icy]=dt_temp*delta_temp;
        }

#define sum_int(var) \
 ( (r_minus*(var[0][0]+var[0][2])+r_plus*(var[2][2]+var[2][0])+ \
   (r_minus*var[0][1]+r_plus*var[2][1]+r_val*(var[1][0]+var[1][2]))+ \
    16.0*r_val*var[1][1])*sqr(h)/24.0 )


      bl_sum=sum_int(bl_int);
      cl_sum=sum_int(cl_int);

      lambda_linear=( (fabs(bl_sum)>almost_same) ? -cl_sum/bl_sum : 0.0 );
      lambda=lambda_linear;  
 
      dchange=dchange-lambda*f_center;

      d7_hold=current_level+time_correct*dchange;
      main_level=old_dd(ix,iy);
      if ((main_level>=0.0)&&(d7_hold<0.0))
       { d7_hold=0.0; }
      else if ((main_level<=0.0)&&(d7_hold>0.0))
       { d7_hold=0.0; }

      d7[ix][iy]=d7_hold;
     }   /* f_center != 0 (delta fn) */

   }    /* looping ix,iy point_ok - mass correction */
  for_all_points
   d4[ix][iy]=d7[ix][iy];   /* update with mass conserved d7 */
}  /* check_mass_conserve */  

/* compute |grad phi| and put in yy[ix][iy] */
void compute_norm_grad(dd,error_width,d_error,max_radius)
double dd[imax][jmax],error_width,*d_error,*max_radius;
{
int ix,iy,k0,k1,k2,k3,k1hold,k2hold,kextra,
    d_count,cnt,i_upwind;
double new_level,current_level,main_level,
       dplus,dminus,dtable,ddtable,dddtable,
       current_radius,bias_level,dchange,out_area,
       dleft,dright,dupwind,
       dxplus,dxminus,dyplus,dyminus,
       temp_error,acmp,bcmp,ccmp;
double second_deriv;

logical slope_limit;

  i_dim=0;
  for_ydomain
   {
    for (ix=4;ix<=mx+2;ix++)
     b00[ix][iy]=(new_dd(ix-2,iy)-new_dd(ix-3,iy))/h;
    for (ix=0;ix<4;ix++)
     b00[ix][iy]=(new_dd1(ix-2,iy)-new_dd1(ix-3,iy))/h;
    for (ix=mx+3;ix<=mx+6;ix++)
     b00[ix][iy]=(new_dd1(ix-2,iy)-new_dd1(ix-3,iy))/h;

    for (ix=0;ix<=mx+5;ix++)
     b01[ix][iy]=(b00[ix+1][iy]-b00[ix][iy])/2.0;
    for (ix=0;ix<=mx+4;ix++)
     b02[ix][iy]=(b01[ix+1][iy]-b01[ix][iy])/3.0;
   }

  i_dim=1;
  for_xdomain
   {
    for (iy=4;iy<=my+2;iy++)
     b10[ix][iy]=(new_dd(ix,iy-2)-new_dd(ix,iy-3))/h;
    for (iy=0;iy<4;iy++)
     b10[ix][iy]=(new_dd1(ix,iy-2)-new_dd1(ix,iy-3))/h;
    for (iy=my+3;iy<=my+6;iy++)
     b10[ix][iy]=(new_dd1(ix,iy-2)-new_dd1(ix,iy-3))/h;

    for (iy=0;iy<=my+5;iy++)
     b11[ix][iy]=(b10[ix][iy+1]-b10[ix][iy])/2.0;
    for (iy=0;iy<=my+4;iy++)
     b12[ix][iy]=(b11[ix][iy+1]-b11[ix][iy])/3.0;
   }

  *max_radius = 1.0;
  *d_error=0.0;
  d_count = 0;

  for_domain
   {
    current_radius=0.0;
    main_level=old_dd(ix,iy);   /* original level */

    /* compute phi_x and phi_y */
    for (i_dim=0;i_dim<2;i_dim++)
     {
      k0=0;
      for (i_upwind=0;i_upwind<=1;i_upwind++)
       {
        if (i_upwind==0)
         {
          k1=k0-1;
          kextra=3;
         }
        else
         {
          k1=k0;
          kextra=1;
         }

        k1hold=-2*k1-1;

        if (i_dim==0)
         dtable=b00[k1+ix+3][iy];
        else
         dtable=b10[ix][k1+iy+3];

        if (i_dim==0)
         { acmp=b01[k1+ix+2][iy];bcmp=b01[k1+ix+3][iy];
           ccmp=b01[kextra+ix][iy]; }
        else
         { acmp=b11[ix][k1+iy+2];bcmp=b11[ix][k1+iy+3];
           ccmp=b11[ix][kextra+iy]; }

        second_deriv=max(fabs(acmp),fabs(bcmp));
        second_deriv=max(second_deriv,fabs(ccmp));
        second_deriv*=2.0;

        if (fabs(acmp)<fabs(bcmp))
         {
          k2=k1-1;k2hold=-k2;ddtable=acmp;
         }
        else
         {
          k2=k1;k2hold=-(k2+2);ddtable=bcmp;
         }

        if (i_dim==0)
         { acmp=b02[k2+ix+2][iy];bcmp=b02[k2+ix+3][iy]; }
        else
         { acmp=b12[ix][k2+iy+2];bcmp=b12[ix][k2+iy+3]; }

        if (fabs(acmp)<fabs(bcmp))
         {
          k3=k2-1;dddtable=acmp;
         }
        else
         {
          k3=k2;dddtable=bcmp;
         }

        dupwind=dtable;
        if (order_space>1)
         dupwind+=ddtable*(k1hold); 
        if (order_space>2)
         {
          /* lack of "h" above taken care of below!! */
          /* allow 3rd order if "radius" > 4 grid cells */
          slope_limit=TRUE;
          if (second_deriv==0.0)
           slope_limit=FALSE;
          else if (fabs(1.0/second_deriv)>4.0)
           slope_limit=FALSE;
          if (! slope_limit)
           dupwind+=dddtable*(k2hold)*(k1hold);
         }

        if (i_upwind==0) dminus=dupwind;
        if (i_upwind==1) dplus=dupwind;
       }   /* looping with i_upwind */
      dright=sign(main_level)*dplus;
      dleft=sign(main_level)*dminus;

      if ((dright<0.0)&&(dleft>0.0))
       current_radius+=max(sqr(dleft),sqr(dright));
      else if (dright<0.0)
       current_radius+=sqr(dright);
      else if (dleft>0.0)
       current_radius+=sqr(dleft);
      else 
       current_radius+=sqr( (dleft+dright)/2.0 );

      if (i_dim==0)
       { dxplus=dplus;dxminus=dminus; }
      else
       { dyplus=dplus;dyminus=dminus; }
     }   /* looping with i_dim */

   current_radius=sqrt(current_radius);

   yy[ix][iy]=current_radius;

   current_level = new_dd(ix,iy);
   if (fabs(current_level)<=error_width)
    {
     *max_radius=max(*max_radius,current_radius);
     d_count++;

     /* normalize error such that 1/2 has same effect as 2 */
     if (current_radius<almost_same)
      temp_error = 1/almost_same;
     else if (current_radius<1.0)
      temp_error = 1/current_radius;
     else
      temp_error = current_radius;

     *d_error+=fabs(1.0-temp_error);
    }   /* current_level close to 0 */
  }   /* looping ix,iy point_ok */

 if (d_count>0) 
  *d_error/=d_count; 

}  /* compute_norm_grad */

void reset_distance1(dd,dd_hld,dnum)
double dd[imax][jmax],dd_hld[imax][jmax];
int dnum;
{
int ix,iy,
    max_steps,steps,cnt,
    step_limit,area_iterate,rk_count;
double time,k,d_error,
       max_radius,new_level,current_level,main_level,
       current_radius,bias_level,dchange,out_area,
       bias,
       time_stop,error_width,
       temp_error,rk_time_reinit, hdivide;
double pn,h0,r_val,r_plus,r_minus;

logical error_small,hld_flag;

  if (! reinit_level)
   return;

  mass_spread=1.0;  /* this insures that the interface will never move */

  for_all_points
  {
   if (! point_ok)
    dd[ix][iy]=gddd(dd,ix,iy);
   d4[ix][iy]=dd[ix][iy]; 
  }
 
  get_current_area(&out_area,dd);
if (use_debug)
  printf("before area,mass_spread,order_space,rk_ord%12.7f %10.5f %5d%5d\n",
    out_area,mass_spread,order_space,rk_max);

  error_width=(density_spread)*h;
 
  time = 0.0;
  time_stop =error_width;
  hdivide=2.0;
  step_limit=(int) (floor((density_spread)*hdivide) );

  if (reinit_test)
   {
    max_steps=(int) (hdivide/h)+1;
    error_width=1.0;
   }
  else if (colortest)
   max_steps=step_limit*40;
  else
   max_steps=-1;

  error_small=FALSE;
  for (steps=0; (! error_small) ;steps++)
    {
     for_all_points  /* this for-loop should not be necessary */
      if (! point_ok)
       d4[ix][iy]=gddd(d4,ix,iy);

     if (max_steps>0)
      if ((steps/10)*10==steps) printf("s%5d\n",steps);

     for (rk_count=0;rk_count<rk_max;rk_count++)
     {

     /* initialize yy[][] with spatial derivative terms */
     compute_norm_grad(dd,error_width,&d_error,&max_radius);

     if ((steps==0)&&(rk_count==0))
      for_domain
       d6[ix][iy]=yy[ix][iy];  /* store for use in mass constraint part */

     /* update in "time"  ; yy hold current value of norm grad phi */
     k=h/hdivide;

     for_domain
       {
        main_level = old_dd(ix,iy);  /* initial level set dd */
        current_level = new_dd(ix,iy);   /* d4 */
        current_radius = grad_dd(ix,iy);  /* yy */
        bias=compute_bias(main_level,&bias_level);

        dchange=bias*(1.0-current_radius);

        if (rk_count==0)
         {
          new_level=current_level+k*dchange;
          d5[ix][iy]=current_level;
          rk_time_reinit=time+k;
         }
        else if (rk_count==1)
         {
          if (rk_max==2)
           {
            new_level=0.5*(d5[ix][iy]+current_level+k*dchange);
            rk_time_reinit=time+k;
           }
          else
           {
            new_level=0.25*(3.0*d5[ix][iy]+current_level+k*dchange);
            rk_time_reinit=time+k/2.0;
           }
         }
        else
         {
          new_level=(2.0*current_level+d5[ix][iy]+2.0*k*dchange)/3.0;
          rk_time_reinit=time+k;
         }

        if ((new_level<0.0) &&
            (main_level>=0.0))
         new_level=0.0;
        else if ((new_level>0.0) &&
                 (main_level<=0.0))
          new_level=0.0;

        d4[ix][iy]=new_level; 
       }    /* looping ix,iy point_ok - initializing the time derivative */

     }  /* looping rk_count */

     time += k; 
     rk_time_reinit=time;

     /* correct d4 to conserve mass - temporary correction in d7 */ 
     if (conserve_mass==TRUE) 
      check_mass_conserve(dd,rk_time_reinit);

     if (max_steps>0)
      {
       if (steps>=max_steps)
        error_small=TRUE;
      }
     else if (time>=time_stop)
      error_small=TRUE;
     else if (steps>=step_limit)
      error_small=TRUE;
     else
      error_small=FALSE;

    }  /* looping steps - until error_small==true */

if (use_debug)
 printf("stp,t,max_t,d_err,max_r,ord%3d%10.5f%10.5f%12.7f%10.5f%5d\n",
     steps,time,time_stop,d_error,max_radius,order_space);
  if (reinit_test)
   {
    steps=0;d_error=0.0;
    for_domain
     {
      max_radius=sqrt(sqr(x[ix])+sqr(y[iy]))-1.0;
      if (fabs(max_radius)<1.0)
       {
        steps++;
        d_error+=fabs(d4[ix][iy]-max_radius);
       }
     }
    printf("distance error of level set cnt,error %5d %10.7f\n",
     steps,d_error/steps);
   }  /* distance error if reinit_test */
  get_current_area(&out_area,d4);
if (use_debug)
  printf("after area%10.5f\n",out_area);

  for_domain
   dd[ix][iy]=d4[ix][iy];
}    /* reset_distance1 */
